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Abstract 

The  development  of  an  algebraic  grid  generation  system  to  track  a solid-liquid  in- 
terface during  directional  solidification  of  a binary  alloy  is  discussed.  A single  mapping, 
constructed  with  tensor  product  B-splines,  is  proposed  for  calculations  of  both  shallow 
and  deep  solidification  cells.  The  initial  spline  coefficients  for  the  coordinate  map- 
ping are  modified  to  minimize  a discrete  functional  that  regulates  the  smoothness  and 
orthogonality  of  the  mesh.  The  use  of  transfinite  blending  function  interpolation  to 
obtain  an  initial  grid  is  examined. 

Keywords:  boundary  fitted  grid  generation,  algebraic  grid  generation,  adaptive,  B- 
splines,  transfinite  blending  functions,  directional  solidification 


1.  Introduction 

One  of  the  well  known  techniques  used  to  study  the  microstructures  that  develop  during 
solidification  of  binary  alloys  is  Bridgman  growth,  a directional  solidification  technique  in 
which  a sample  of  the  alloy  is  drawn  through  a constant  temperature  gradient  at  a uniform 
rate  of  speed,  V,  as  shown  in  Figure  1.  A considerable  amount  of  theoretical  work  has 
focused  on  examining  the  morphological  instabilities  of  the  growing  soHd-liquid  interface  [1- 
6].  MulHns  and  Sekerka  [1]  used  Hnear  stability  theory  to  predict  the  critical  velocity  for  the 
onset  of  instabihty  for  a planar  interface.  Experimental  observations  confirm  the  validity  of 
their  results  and  show  that  after  the  onset  of  instability  the  structure  of  the  interface  can 
change  from  planar  to  cellular  to  dendritic  and  back  again  as  the  growth  velocity  is  increased 
[7,  8].  Coriell  et  al.  [2]  extended  the  results  of  Mulhns  and  Sekerka,  including  the  effects 
from  convection  in  the  Hquid. 

Although  cellular  microstructures  are  much  simpler  than  dendritic,  as  the  control  param- 
eters, growth  velocity  or  temperature  gradient,  are  changed,  the  cells  may  become  very  deep 
and  naxrow  with  re-entrant  bulb-Hke  shapes.  To  successfully  track  the  interface,  the  grid 
generation  mapping  must  adapt  to  large  deformations  of  the  interface  shape  while  maintain- 
ing as  much  orthogonality  and  smoothness  as  possible.  Ettouney  and  Brown  [4]  successfully 
modeled  slightly  nonplanaj  interfaces  by  using  cin  algebraic  grid  generation  system  where  the 
interface  was  described  in  terms  of  a univariate  function.  Unfortunately,  this  representation 
of  the  interface,  called  a Monge  transformation,  fails  for  cells  with  vertical  or  re-entrant 
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Figure  1:  Bridgman  growth  technique. 


walls.  To  overcome  that  problem,  Ungar  and  Brown  [5]  developed  a representation  defined 
by  the  division  of  the  interface  into  disjoint  sections  that  could  be  expressed  as  separate 
Monge  transformations  written  in  either  polar  or  rectangular  coordinates.  With  this  mixed 
transformation,  Ungar  and  Brown  modeled  cellular  interfaces  with  grooves  as  much  as  15 
times  longer  than  the  cell  wavelength.  Ideally,  analysts  want  a single  grid  generation  tech- 
nique that  can  track  the  interface  as  it  changes  from  a shallow  deformation  to  a narrow  and 
deeply  grooved  cell  with  re-entrant  sidewalls.  Furthermore,  it  should  capture  the  change  in 
wavelength  if  the  cell  sphts. 

Tsiveriotis  and  Brown  [6]  have  made  some  progress  in  this  area  by  using  a two-step 
procedure.  In  one  direction  the  coordinate  is  defined  by  a generalized  Poisson  equation  with 
a scahng  condition  that  forces  the  grid  lines  into  concave  areas  of  the  interface.  The  other 
coordinate  is  obtained  by  minimizing  smoothness  and  orthogonality  functionals  similar  to 
those  of  Brackbill  and  Saltzman  [9].  Tsiveriotis  and  Brown  show  impressive  examples  of 
grids  developed  for  narrow  and  re-entrant  cells,  but  it  is  not  clear  whether  this  method  is 
as  successful  as  the  previous  method  of  Ungar  and  Brown  because  the  most  narrow  cells  in 
the  earher  paper  are  not  redone.  Also,  to  obtain  better  results,  they  fix  a horizontal  line 
near  the  interface  to  decouple  the  far  field  domain  from  the  domain  near  the  interface.  This 
essentially  creates  two  separate  domains  whose  grid  cells  are  not  smoothly  connected  near 
the  line.  This  is  fine  for  finite  element  calculations,  but  suggests  problems  for  those  who  are 
interested  in  using  the  technique  to  create  grids  suitable  for  finite  difference  schemes. 

This  paper  describes  progress  to  date  in  the  development  of  an  algebraic  technique  for 
generating  a boundary/interface  fitted  coordinate  system  that  is  smooth  enough  to  be  used 
for  finite  difference  calculations.  No  partial  differential  equations  (pdes)  are  solved  to  obtain 
the  coordinate  system,  and  the  same  system  is  used  for  the  entire  domain.  The  computation 
of  the  interface  is  not  discussed  in  this  paper.  It  is  assumed  that  the  interface  is  available 
as  a discrete  set  of  points,  having  been  determined  by  some  means  such  as  an  evaluation  of 
the  Gibbs-Thomson  equation  which  relates  the  interface  mean  curvature,  concentration  of 
solute  in  the  liquid,  and  the  melting  temperature  near  the  interface  [8,  10].  The  feasibility 
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of  the  grid  generation  technique  is  evaluated  by  examining  its  abiUty  to  create  grids  that 
conform  to  interface  shapes  typical  of  those  seen  in  experimental  observations  of  directional 
sohdification  of  binary  alloys  by  the  Bridgman  growth  technique  [7,  8].  The  grid  generation 
algorithm  is  an  extension  of  an  algebraic  technique  for  generating  boundary  fitted  grids  [11]. 
This  paper  discusses  the  modifications  needed  and  addresses  the  problems  in  tracking  an 
interface  that  deforms  from  a planar  shape  to  a deep,  re-entrant  cellular  microstructure. 

Section  2 provides  a brief  description  of  two  algebraic  techniques  for  generating  coordinate 
systems:  transfinite  blending  function  interpolation,  the  technique  used  to  obtain  an  initial 
approximation  for  the  interface  fitted  mesh,  and  a boundary  fitted  grid  generation  technique 
that  uses  a mapping  constructed  with  B-splines.  The  B-spHne  mapping  is  extended  to 
create  the  mapping  for  the  interface  tracking  grid  generation  system  . Section  3 explains  the 
construction  of  the  grid  generation  mapping  and  Section  4 examines  the  results  to  date. 


2.  Algebraic  Grid  Generation  Techniques 

In  algebraic  grid  generation,  a direct  transformation  describes  the  relationship  between  the 
computational  and  physical  domains.  The  transformation  is  constructed  so  that  it  interpo- 
lates the  boundary  points  and/or  points  in  the  interior.  It  may  also  be  constrained  to  match 
derivatives.  No  partial  differential  equations  are  solved  to  obtain  the  curvilinear  coordinates, 
so  algebraic  techniques  can  be  easier  to  construct  than  pde  methods,  and  give  easier  con- 
trol over  grid  characteristics  such  as  orthogonality  and  grid  point  spacing.  However,  these 
methods  are  sometimes  criticized  for  allowing  discontinuities  on  the  boundary  to  propagate 
into  the  interior  and  for  not  generating  grids  as  smooth  as  those  generated  by  pde  methods. 
Nevertheless,  algebraic  techniques  have  been  used  successfully  to  generate  grids  in  both  two 
and  three  dimensions  [11-16]. 


2.1.  Transfinite  Blending  Function  Interpolation 

One  of  the  most  common  and  easiest  algebraic  methods  to  implement  is  transfinite  blending 
function  interpolation,  where  interior  regions  are  represented  in  terms  of  boundary  functions. 
This  technique,  which  has  been  widely  used  for  problems  in  grid  generation  [11, 12, 14, 15, 16], 
was  originally  developed  for  problems  in  computer  aided  design.  Gordon  and  Hall  [14] 
adapted  and  apphed  the  technique  to  grid  generation  for  finite  element  and  finite  difference 
calculations  in  the  late  sixties  and  early  seventies.  A simple  example  is  illustrated  by  the 
mapping  T from  the  unit  square  I2  to  the  physical  domain  defined  by 


Tu,r/)  = 

i 

3 

- (1) 


where  0 = 6 < •••  < ^Af  = 1 and  0 = rjo  < . . . < Tjpf  = 1.  The  continuous  vector- 
valued function  f describes  the  boundary  of  the  physical  domain.  $ and  ^ are  “blending” 
or  “connecting”  functions  that  satisfy  the  conditions 
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1,  i = j 
0>  j 


(2) 


^kim)  = { 5’  (3) 

Therefore,  if  M=N=1,  one  might  choose  $ and  ^ to  be  linear  Lagrange  polynomials  so  that 
T matches  f on  the  boundary  of  l2-  If  M and  N are  larger  and  derivative  information  is 
given,  the  blending  functions  can  be  chosen  so  that  interior  curves  and  the  derivatives  there 
are  matched.  Unfortunately,  transfinite  blending  function  interpolation  may  allow  boundary 
singularities  to  propagate  into  the  interior  of  the  mesh.  Furthermore,  if  the  boundary  is 
nonconvex,  the  grid  lines  may  overlap  the  boundary.  The  choice  of  blending  functions  and 
parametrization  of  the  boundary  function  are  very  important  in  alleviating  these  problems. 
Transfinite  blending  function  interpolation  provides  a relatively  easy  way  of  obtaining  an 
initial  grid  that  can  be  refined  and  smoothed  by  other  techniques,  whether  algebraic,  pde, 
or  variational.  However,  unfortunately,  some  techniques  cannot  use  the  grid  if  it  contains 
areas  where  the  grid  lines  overlap,  that  is,  areas  of  negative  Jacobian  [17]. 


2.2.  Boundary  Fitted  Grid  Generation  Using  Tensor  Product  B-splines 


Saunders  [11]  developed  an  algebraic  grid  generation  system  that  uses  a mapping  T from 
the  unit  square  I2  to  a physical  domain  of  arbitrary  shape  defined  by 


I y{i,v)  j \ Er=i  ) ’ 


(4) 


where  0 and  each  Bij  is  the  tensor  product  of  cubic  B-sphnes.  Hence,  Bij^^^r})  = 

Bi{i)Bj{r))  where  Bi  and  Bj  are  elements  of  cubic  B-spline  sequences  associated  with  finite 
nondecreasing  knot  sequences,  say,  and  respectively.  As  defined  by  de 

Boor  [18],  B-sphnes  are  essentially  piecewise  polynomials  with  continuity  conditions  at  each 
breakpoint  determined  by  the  repetition  of  the  breakpoint  value  in  the  associated  knot 
sequence.  For  a typical  cubic  B-spline  Bj,  the  value  of  the  B-spHne  is  determined  by  the  five 
knots  tj,  tj+i,  tj+2,  tj+3,  tj+4.  Its  support  is  small,  that  is,  Bj  can  be  nonzero  only  on  the 
interval  [tj,tj+4\.  Consequently,  only  four  B-splines,  Bj-3,  Bj  can  be  nonzero  on 

the  interval  (tj,tj^i).  The  initial  coefficients  are  chosen  so  that  the  mapping  approximates 
transfinite  blending  function  interpolation.  More  specifically,  the  coefficients  are  selected  to 
produce  a variation  diminishing  spHne  approximation  to  the  transfinite  blending  function 
interpolant  T defined  in  equation  (1).  In  short,  this  means  the  coefficients  are  obtained  by 
evaluating  T at  average  knot  values  as  discussed  in  [18].  This  shape  preserving  approximation 
reproduces  straight  lines  and  preserves  convexity  [11,  18]  . To  increase  the  orthogonality  and 
smoothness  of  the  grid  lines  the  initial  coefficients  are  modified  to  minimize  the  functional 


where  T denotes  the  grid  generation  mapping,  J is  the  Jacobian  of  the  mapping,  and  wi  and 
W2  are  weight  constants.  When  is  large,  the  variation  of  the  Jacobian  values  at  nearby 
points  will  be  small,  thereby  decreasing  skewness.  When  W2  is  large,  the  dot  product  term 
will  be  small,  causing  the  grid  fines  to  approach  orthogonality.  To  avoid  solving  the  Euler 
equations  for  the  variational  problem,  this  functional  is  approximated  in  the  computer  code 
by  the  sum 


Ji+i,j  Jjj 

H . 


+ 


Af  At; 
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(6) 


+ Y^W2Dot^jAiArj 

where  Jij  is  the  Jacobian  value  and  Dotij  is  the  dot  product  of  dTl  j and  at  mesh 

point  {(nVj)  square.  G is  actually  a fourth  degree  polynomial  in  each  spline 

coefficient  so  the  minimum  is  found  by  using  a cyclic  coordinate  descent  technique  which 
sequentially  finds  the  minimum  with  respect  to  each  coefficient.  This  technique  allows  the 
minimization  routine  to  take  advantage  of  the  small  support  of  B-sphnes  when  evaluating 
the  sums  that  comprise  G. 

An  application  of  the  grid  generation  algorithm  is  shown  in  Figures  2 and  3 for  a puzzle 
shaped  domain.  The  boundary  and  initial  grid  are  shown  in  Figure  2.  The  initial  grid 
was  constructed  using  linear  Lagrange  polynomials  for  the  blending  functions.  Note  that 
the  grid  lines  overlap  the  nonconvex  boundary.  Figure  3a  shows  the  mesh  obtained  after 
the  spline  coefficients  are  modified  to  minimize  the  smoothing  functional.  The  skewed  areas 
have  been  eliminated  and  the  overlapping  grid  fines  have  been  pulled  into  the  interior.  The 
refined  mesh  in  Figure  3b  is  obtained  by  evaluating  T at  additional  points  on  the  square. 
The  concentration  near  the  bottom  is  achieved  by  applying  an  exponential  function  to  the  rj 
variable.  The  algorithm  discussed  in  this  paper  extends  the  boundary  fitted  mapping  to  fit  a 
curve  in  the  interior  of  the  physical  domain.  As  in  the  case  of  the  boundary  fitted  algorithm, 
it  will  be  shown  that  the  new  algorithm  is  robust  enough  to  smooth  and  untangle  an  initial 
grid  containing  skewed  and  overlapping  grid  fines. 


(a)  Puzzle  shaped  boundary.  (b)  Initial  grid. 

Figure  2:  Puzzle  shaped  domain.  Initial  grid  was  produced  using  an  approximation  of 
transfinite  blending  function  interpolation. 


3.  An  Interface  Tracking  Grid  Generation  System 

The  boundary  fitted  grid  generation  mapping  discussed  in  the  previous  section  forms  the 
basis  for  the  interface  tracking  mapping.  However,  the  mapping  must  now  match  the  interface 
curve  on  the  interior  of  the  physiccil  domain  in  addition  to  fitting  the  outer  physical  boundary. 
Furthermore,  the  system  must  be  adaptive  since  the  grid  fines  must  change  to  follow  the 
deforming  interface  while  maintaining  as  much  smoothness  and  orthogonality  as  possible. 
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(a)  Optimized  grid. 


(b)  Optimized  grid  refined. 


Figure  3:  Optimized  puzzle  grids.  Minimization  of  smoothing  functional  successfully  pulls 
grid  lines  inside  the  boundary. 


3.1.  Grid  Generation  Mapping 

The  proposed  grid  generation  mapping,  T,  maps  the  unit  square,  J2)  onto  the  physical 
domain  and  is  constructed  so  that  the  interface  is  the  coordinate  curve  77  = 1/2  as  shown  in 
the  figure.  As  before,  the  mapping  has  the  form 

~ I y(^. j “ I Er=i  E"=i PML ri))' 

where  0 < ^,77  < 1 and  Bij{(yr))  = Bi(()Bj{T})  where  Bi  and  Bj  are  elements  of  cubic  B- 
spline  sequences  associated  with  finite  nondecreasing  knot  sequences,  and 

respectively.  The  spline  coefficients  for  T can  be  divided  into  three  groups.  The  boundary 
coefficients  are  the  coefficients  of  those  Bij  that  are  nonzero  on  the  boundary  of  I2.  Since  T 
is  defined  so  that  the  interface  corresponds  to  coordinate  curve  77  = 1/2,  the  coefficients  of 
the  Bij  that  are  nonzero  when  77  = 1/2  are  called  the  interface  coefficients.  The  remaining 
coefficients  are  called  the  interior  coefficients. 

Initially,  the  coefficients  are  chosen  to  approximate  the  following  transfinite  blending 
function  interpolant  that  matches  the  outer  boundary  and  interface  of  the  physical  domain: 


t=0 


j=0 

+ 1/2) 

t“0  j=0 

^ df 

- E^.(0/3(’?)^te.i/2) 


(8) 


where  = 0,  = 1 and  770  = 0,  771  = 1/2,  772  = 1.  The  continuous  vector  valued  function 

f,  constructed  from  boundary  data  input  by  the  user,  maps  the  boundary  of  the  square  to 
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Figure  4:  Grid  generation  mapping. 

the  boundary  of  the  physical  domain  and  the  line  77  = 1/2  to  the  interface.  The  transfinite 
mapping  T interpolates  f on  the  boundary  and  matches  f and  dfldr)  at  rj  = 1/2.  The 
blending  functions  are  linear  Lagrange  polynomials  satisfying 

$i(6)  = for  = o>i 

while  and  (3  are  hermite  quintic  blending  functions  satisfying 


= Srs 

'^r{Va)  = 0 for  r,s  = 0, 1, 2 


PiVs)  = 0 

/3X'n»)  = for  s = 0,1,2  . 


Although  linear,  quadratic,  and  cubic  polynomials  were  also  tested  as  blending  functions  for 
the  T)  coordinate,  the  hermite  quintic  polynomials  produced  the  least  amount  of  skewness 
and  overlap  of  grid  cells.  To  force  orthogonality  at  the  interface,  the  components  of  the 
derivative 


drj 


(^.1/2) 


are  chosen  to  be 


where 


L = 


(dh\\(^y 

and  A"  is  a user  defined  orthogonality  constant  that  regulates  the  magnitude  of  the  normal 
vectors  at  the  interface.  The  initial  coefficients  for  the  grid  generation  mapping  T are  defined 
by 


a 


= for  t = j = 1, 


,n 


(9) 


where  5*  = (-St+i  + . . . + 3j+3)/3,  i = 1, . . . , m and  + . . . + j = 1, . . . , n. 

With  these  coefficients  T produces  a variation  diminishing  spline  approximation  to  T.  Vari- 
ation diminishing  spHne  approximations,  discussed  in  detail  in  [18],  are  shape  preserving 
approximations  that  reproduce  straight  lines  and  preserve  convexity. 

As  with  the  boundary  fitted  system,  the  smoothness  and  orthogonality  of  grid  Hnes  is 
enhanced  by  modifying  the  coefficients  to  minimize  the  discrete  smoothing  functional  de- 
fined in  equation  (6).  However,  the  flexibihty  of  the  boundary  and  the  interface  coefficients 
is  hmited  since  even  small  changes  in  the  coefficients  can  destroy  the  shape  preserving  prop- 
erties of  the  spline  mapping.  Fortunately,  the  B-spline  knot  sequences  {5}  and  {t}  can  be 
chosen  so  that  the  boundary  coefficients  affect  only  the  boundary  ajid  a narrow  area  inside 
the  boundary.  Using  the  continuity  properties  of  B-spHnes  with  repeated  knots  [18],  one  can 
show  that  if  the  first  four  knots  of  and  are  0,  the  last  four  knots  are  1,  and 

the  rest  located  in  the  interval  (0,1),  then  the  only  boundary  coefficients  are  aij,/3ij  and 
GLmjiPmj  j = 1 , . . . , 71  and  and  OLim^in  for  i = 1, . . . ,r7i.  Unless  all  the  “interior” 

knots  are  clustered  near  the  center  of  (0, 1),  the  boundary  coefficients  will  have  a minimal 
effect  on  the  interior  of  the  square.  For  example,  aij  and  Pij  will  only  affect  points  on  the 
support  of  Bijj  that  is,  the  narrow  band  [0,35]  x illustrated  in  Figure  5,  where 

S5  is  the  first  interior  knot  in  In  fact  the  area  of  influence  of  all  the  boundary 


Figure  5:  Support  of  tensor  product  B-spline  Bij. 
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coefficients  will  be  a narrow  region  along  the  boundary  of  the  square.  Consequently,  even  if 
the  boundary  coefficients  are  fixed  to  guarantee  that  “boundary  fittedness”  is  maintained, 
the  smoothing  functional  should  be  able  to  produce  a significant  amount  of  orthogonahty 
and  smoothness  in  the  grid. 

Ideally,  however,  one  would  like  to  include  the  boundary  coefficients  in  the  minimization 
process  to  ensure  that  grid  lines  near  the  boundary  are  smooth  and  orthogonal.  When 
the  boundary  is  rectangular,  as  it  is  with  Bridgman  growth,  this  is  indeed  possible.  The 
four  sides  are  vertical  or  horizontal  line  segments  satisfying  either  x{0,t])  = xq,  x[l,rj)  = xi, 

0)  = 2/0,  or  ?/(f , 1)  = ^1,  where  xq,  xi,  yo,  and  yi  are  constants.  By  fixing  the  appropriate 
a oi  (5  coefficients  one  can  guarantee  that  these  boundary  equations  will  remain  satisfied. 
For  example,  on  the  left  vertical  boundary  where  cc(0,7/)  = xq,  only  Bij,j  = l,...,n  are 
nonzero  there  so  that  x{0^r})  = Yij=i  Hence,  ^(O,?/)  = Xq  remains  true  if  aij, 

j = 1, . . . ,n,  that  is,  the  a boundary  coefficients  for  that  side,  are  fixed.  The  corresponding 
P coefficients,  /9ij,  j = 1, . . . , n,  are  free  to  move.  Similarly,  for  the  right  side,  the  a boundary 
coefficients,  amj,  j = 1, . . . , n are  fixed.  For  the  horizontal  sides,  the  /3  boundary  coefficients 
are  fixed.  With  the  boundary  coefficients  fixed  in  this  manner,  minimizing  the  smoothing 
functional  causes  a redistribution  of  the  boundary  points  along  the  vertical  or  horizontal 
lines.  Of  course,  it  is  possible  that  a boundaxy  point  could  be  pushed  beyond  the  endpoint 
of  the  line  segment,  but  this  produces  a negative  Jacobian,  a problem  generally  corrected  by 
the  smoothing  action  of  the  Jacobian  terms  in  the  functional. 

Since  the  interface  may  be  quite  complex,  very  Httle  change  is  permitted  in  the  interface 
coefficients  when  the  other  coefficients  are  adjusted  by  the  minimization  of  the  smoothing 
functional.  Unfortunately,  unlike  the  boundary  coefficients,  the  interface  coefficients  affect 
a significant  number  of  points  on  the  interior  besides  those  mapped  onto  the  interface.  To 
determine  the  interface  coefficients  one  first  finds  I such  that  1/2  E Then  the 

interface  coefficients  will  be  l3ik  where  1 < i < m and  I — 3 < k < 1.  These  are  the 
coefficients  of  the  tensor  product  B-splines  that  might  be  nonzero  when  77  = 1/2.  However, 
fixing  these  coefficients  affects  the  mapping  T not  only  on  [0, 1]  x but  also  on  the 

much  larger  bajid,  [0,1]  x [ti_3,tf+4]  seen  in  Figure  6,  that  is,  the  total  support  of  all 
the  tensor  product  B-splines  that  are  possibly  nonzero  on  [0, 1]  x This  band  can  be 

narrowed  significantly  by  using  a larger  number  of  knots  for  the  t sequence  and  concentrating 
some  near  1/2.  However,  concentrating  the  knots  too  closely  will  affect  the  smoothness  of 
the  grid  fines.  On  the  other  hand,  since  a continuous  B-spfine  is  close  to  zero  near  the 
boundary  of  its  support,  all  interface  coefficients  do  not  equally  affect  the  mapping  of  the 
interface.  In  particular,  allowing  the  coefficients  and  for  1 < i < m to  move 

appears  to  have  very  little  effect  on  the  accuracy  of  the  interface  mapping.  Nevertheless,  the 
inflexibility  of  the  interface  coefficients  suggests  that  they  should  be  chosen  to  produce  as 
much  orthogonality  and  smoothness  as  possible  at  the  outset. 

3.2.  The  Algorithm 

Although  the  calculation  of  the  interface  must  be  addressed  in  order  to  integrate  the  grid 
generation  system  with  the  equations  that  model  the  directional  solidification  problem,  the 
current  code  assumes  that  both  the  boundary  and  interface  data  are  available  as  sets  of 
discrete  points.  The  focus  here  is  on  demonstrating  the  ability  of  the  grid  generation  mapping 
to  adapt  as  the  interface  changes. 

Using  the  boundary  data  and  initial  interface  data,  the  code  first  computes  the  B-splines 
needed  to  define  T.  K and  are  the  sequences  chosen  to  define  the  B-splines, 
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Figure  6:  Support  of  tensor  product  B-splines  nonzero  when  r/  = 1/2. 

then  the  total  number  of  coefficients  used  in  the  mapping  will  be  2 x m x n.  Consequently, 
although  the  number  of  knots  must  be  sufficient  to  accurately  capture  the  shape  of  a com- 
plicated interface,  this  must  be  balanced  by  the  fact  that  increasing  the  number  of  knots 
increases  the  size  of  the  optimization  problem.  Fortunately,  the  accuracy  of  the  mapping 
is  also  improved  by  concentrating  the  t knots  near  the  interface.  This  is  accomplished  by 
replacing  {t}  with  the  sequence  {7(t)}  where  7 is  the  hyperbolic  sine  function 

sinh(2ct  — c)  -|-  sinh(c) 

""  2sinh(c) 

where  c is  a constant  that  determines  the  degree  of  concentration.  This  has  the  additional 
beneficial  effect  of  decreasing  the  area  influenced  by  the  interface  coefficients.  The  mesh  over 
which  the  smoothing  functional  is  evaluated  must  have  at  least  2 x m x n grid  points,  and 
the  concentration  of  points  in  the  rj  direction  should  be  similar  to  the  concentration  of  knots 
for  sequence  {t}. 

As  stated  above,  the  transfinite  blending  function  mapping  T is  evaluated  at  average 
knot  values  to  obtain  coefficients  that  produce  a variation  diminishing  spline  approximation 
to  T.  The  coefficients  are  modified  to  minimize  the  discrete  smoothing  functional  to  obtain 
more  smoothness  and  orthogonality  in  the  grids.  This  step  is  unnecessary  if  the  interface 
is  planar.  The  actual  grid  is  obtained  by  evaluating  T at  equally  spaced  values  of  f and 
Tj.  Replacing  7 with  the  sinh  mapping  7(77)  concentrates  the  grid  near  the  interface.  To 
create  a grid  for  the  next  interface,  the  code  first  reads  the  data  for  the  new  interface.  The 
old  interface  coefficients  plus  one  layer  of  coefficients  above  and  below  them,  that  is, 

Pik  where  1 < i < m and  / — 4 < k < / +!,  are  replaced  by  coefficients  that  produce 
a variation  diminishing  spline  approximation  to  the  transfinite  blending  function  mapping 
that  interpolates  the  new  interface.  The  smoothing  functional  is  then  minimized  in  order  to 
eliminate  any  overlap  of  grid  lines  and  improve  smoothness  and  orthogonahty.  The  process 
is  repeated  to  obtain  a grid  for  the  next  interface. 

These  steps  are  illustrated  in  Figure  7 which  shows  a sequence  of  grids  concentrated 
neaj  curves  that  demonstrate  the  type  of  sinusoidal  deformation  that  appears  soon  after  the 
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onset  of  instability  during  the  directional  solidification  of  a metal  alloy.  The  first  grid,  rep- 
resenting a planar  interface,  requires  no  additional  smoothing.  The  eight  other  grids  consist 
of  pairs  representing  the  initial  and  smoothed  grids,  respectively.  The  minimization  routine 
successfully  “untangles”  overlapping  Hnes  while  enhancing  the  smoothness  and  orthogonality 
of  the  grid  lines.  A concentration  of  grid  points  near  the  interface  was  obtained  by  evaluating 
with  sinh  concentration  constant  c = 4. 

3.3.  Implementation 

The  usefulness  of  this  technique  for  generating  grids  for  solidification  problems  will  largely 
be  determined  by  how  fast  and  efficiently  the  grids  can  be  produced.  Of  course,  the  primary 
bottleneck  in  the  grid  generation  routine  is  the  minimization  of  the  smoothing  functional.  De- 
pending on  the  number  of  knots  used  to  define  the  grid  generation  mapping,  the  optimization 
problem  can  be  of  fairly  large  scale.  For  example,  if  the  knot  sequences  {3}  and  {t}  both  have 
22  elements,  then  the  number  of  spHne  coefficients,  aij  and  ^ij,  will  be  2 x (22  — 4)  x (22  — 4) 
or  648  [18].  Hence  we  have  a minimization  problem  involving  648  variables.  An  ongoing 
problem  is  to  determine  the  best  optimization  technique  to  use.  Currently,  the  code  uses  a 
modified  cycHc  coordinate  descent  algorithm  which  minimizes  a multivariable  function  by  se- 
quentially finding  the  minimum  with  respect  to  each  variable.  The  cyclic  coordinate  descent 
method  is  relatively  easy  to  implement  and  requires  no  gradient  information  to  determine 
the  direction  of  descent,  but  in  general  it  has  a much  slower  convergence  rate  than  conjugate 
gradient  or  quasi-Newton  methods  [19].  However,  the  simplicity  of  the  technique  makes  it 
easy  to  exploit  the  properties  of  B-splines  during  the  minimization  process.  For  each  spline 
coefficient,  the  sum  in  equation  (6)  need  only  be  computed  over  the  cirea  that  the  coefficient 
affects.  The  small  support  of  B-splines  means  that  this  area  can  be  quite  small.  The  example 
illustrated  in  Figure  8 shows  the  support  of  tensor  product  B-spline  Bej  for  specified  knot 
sequences  {s}  and  {t}  and  a given  distribution  of  mesh  points  on  the  computational  domain. 
If  the  functional  is  minimized  with  respect  to  coefficients  asj  or  I3ej  then  the  sums  in  (6) 
need  only  be  computed  over  i,  j such  that  2 < z < 5 and  4 < j < 7.  The  B-spline  values  and 
their  derivatives  are  computed  using  the  de  Boor  routines  BVALUE,  BSPLVN,  BSPLVD, 
and  INTERV  [18]. 

Even  with  the  smaller  sums  a significant  amount  of  time  can  be  wasted  during  the 
execution  of  the  minimization  routine  if  the  code  can  not  take  advantage  of  the  optimization 
features  of  the  computer  being  used.  At  the  present  time  the  code  is  running  on  a Cray  Y- 
MP4E/232  computer  system  having  two  central  processing  units  (CPUs).  The  small  support 
of  tensor  product  B-splines  facilitates  the  design  of  an  algorithm  that  takes  advantage  of  the 
vectorization  features  of  the  computer.  One  can  easily  show  that  the  support  of  a tensor 
product  B-spline  Bij  will  be  disjoint,  except  possibly  on  the  support  boundary,  from  the 
support  of  any  tensor  product  B-spline  Bpq  where  either  |p  — fl  > 4 or  \q  — j\  >4.  This  is 
illustrated  in  Figure  9 for  Bsj.  This  means  that  the  minimization  code  can  be  designed 
so  that  the  smoothing  functional  is  minimized  with  respect  to  the  coefficients  associated 
with  those  tensor  product  B-splines  having  disjoint  support  simultaneously.  The  sums  must 
be  restricted  so  that  they  are  computed  only  over  those  points  that  faU  inside  the  support 
area,  but  this  does  not  appear  to  decrease  the  effectiveness  of  the  minimization  routine. 
Rewriting  the  code  to  minimize  a vector  of  coefficients  simultaneously  produced  a speedup 
by  a factor  close  to  4.5.  By  integrating  small,  frequently  called  routines  into  the  calling 
routines  (inlining)  and  eliminating  redundant  code,  a larger  speedup  was  obtained.  For 
example,  the  original  scalar  code  eliminated  the  overlap  of  grid  lines,  that  is,  areas  with 
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Figure  7:  Grids  simulating  the  gradueil  deformation  of  an  interface. 
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Figure  8:  Support  of  tensor  product  B-spline  Bqj. 
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Figure  9:  Tensor  product  B-splines  with  disjoint  support.  The  shaded  area  shows  the  support 
of  Bsj.  The  surrounding  boxes  show  other  tensor  product  B-splines  with  disjoint  support. 
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Table  1:  User  CPU  Time  (seconds) 


Iterations 

Scalar 

Vectorized 

Vectorized-f 

5 

216.9 

48.6 

13.3 

10 

442.36 

97.04 

27.05 

negative  Jacobian,  by  determining  the  appropriate  interval  that  would  guarantee  a positive 
Jacobian  for  each  coefficient.  However,  this  is  unnecessary  since  the  Jacobian  terms  of  the 
smoothing  functional  tend  to  perform  the  same  function  by  untanghng  grid  lines  to  minimize 
the  difference  in  Jacobians  at  nearby  grid  points. 

Table  1 shows  the  user  CPU  times  for  a problem  of  moderate  size.  For  this  problem,  m = 
18  and  n ~ 24  so  that  the  number  of  spline  coefficients  is  2 x 18  x 24  = 864.  The  minimization 
was  performed  on  a 60  x 50  mesh.  Since  the  system  contains  only  two  processors,  very  little 
effort  was  focused  on  paxallelizing  the  code.  The  automatic  parallelization  option  of  the 
compiler  was  tested,  but  it  produced  only  a minor  speedup  in  the  code.  The  table  shows  the 
single  processor  times  for  five  and  ten  iterations  of  the  minimization  routine  for  three  cases: 
the  scalar  code,  the  initial  vectorized  version  where  the  code  was  rewritten  to  minimize  with 
respect  to  a vector  of  coefficients  (vectorized),  and  the  improved  vectorized  version  in  which 
some  routines  were  inlined  and  others  eliminated  (vectorized-f ).  The  final  vectorized  run  is 
over  16  times  faster  than  the  scalar  run. 


4.  Results  and  Discussion 

For  all  the  examples  shown,  the  number  of  spline  coefficients  used  to  define  the  grid  gen- 
eration mapping  is  2 x m x n where  m is  the  number  of  B-spHnes  associated  with  knot 
sequence  {s},  or  the  “f”  sequence,  and  n is  the  number  associated  with  knot  sequence  {t}, 
the  sequence.  Recall  that  a concentration  of  grid  points  near  the  interface  (coordinate 
curve  7)  = 1/2)  is  obtained  by  evaluating  T at  (^,7(77)).  To  increase  the  accuracy  of  the 
spline  mapping  and  decrease  the  area  affected  by  the  interface  coefficients,  the  sinh  function 
is  also  used  to  concentrate  the  77  knots  near  the  interface.  This  means  that  the  sequence 
{7(0}  used  instead  of  {t}  in  the  definition  of  the  mapping.  The  constant  values  of  c 
associated  with  these  two  t5rpes  of  concentration  are  called  the  mesh  concentration  and  the 
knot  concentration,  respectively.  The  orthogonality  constant  K represents  the  magnitude  of 
the  normal  vectors  at  the  interface  as  defined  in  section  3.1.  It  should  also  be  noted  that  the 
term  “initial  grid”  is  used  to  refer  to  any  grid  computed  using  the  initial  spline  coefficients. 

In  the  first  example  the  bottom  of  the  puzzle  boundary  displayed  in  Figure  2 is  now 
shown  as  an  interface.  This  is  similaj  to  the  re-entrant  shape  that  commonly  appears  in 
cellular  microstructures.  For  this  example  m = 19  and  n = 20.  The  orthogonality  constant 
K is  ^ and  the  knot  concentration  is  3.  The  variation  diminishing  spline  approximation  to 
T produces  the  mesh  shown  on  the  left  in  Figure  10.  Although  the  grid  cells  axe  skewed  in 
some  axeas,  the  grid  appears  orthogonal  near  the  interface.  Hence,  the  initial  grid  already 
looks  fairly  good  there.  This  is  important  since  most  of  the  interface  coefficients  remain  fixed 
throughout  the  minimization  process.  Using  a 40  x 42  initial  grid  with  mesh  concentration 
3,  the  optimization  routine  significantly  improves  the  smoothness  and  general  orthogonality 
of  the  grid  cells.  This  is  illustrated  in  the  grid  on  the  right,  computed  after  forty  iterations. 
The  Jacobian  and  orthogonality  weights  used  for  the  smoothing  functional,  W\  cind  ii;2,  were 
0.65  and  10,  respectively. 
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(a)  Initieil  grid. 


(b)  Optimized  grid. 


Figure  10:  Grids  produced  using  linear  and  hermite  quintic  blending  functions.  Grid  on 
right  was  produced  after  forty  iterations  of  the  smoothing  routine. 

To  illustrate  the  flexibility  of  the  grid  generation  routine,  the  next  six  grids  show  grids 
that  fit  curves  that  are  typical  of  the  shapes  seen  during  deformation  of  an  interface  from  a 
sinusoidal  shape  to  a deep  cell.  These  grids  actually  continue  the  sequence  that  was  started 
in  Figure  7.  The  first  grid  in  Figure  11  is  the  final  grid  in  Figure  7.  The  grids  show 
the  meshes  obtained  after  the  smoothing  functional  is  minimized.  The  initial  grids,  most 
of  which  contciined  overlapping  grid  fines,  are  not  shown.  For  all  of  the  grids  m = n = 18 
and  the  knot  concentration  is  3.  The  minimization  was  performed  on  a 40  x 40  mesh  using 
a mesh  concentration  of  3.  An  orthogonality  constant  of  AT  = 4 was  used  for  the  grids  in 
Figure  11  and  Figure  12a.  This  was  increased  to  4.5  for  Figure  12b  and  to  5 for  the  grids  in 
Figure  13  to  maintain  the  smoothness  and  orthogonality  near  the  interface  as  the  interface 
cell  deepened.  In  every  case  the  smoothing  routine  successfully  removed  overlapping  grid 
fines.  The  number  of  iterations  of  the  minimization  routine  executed  varied  from  5 or  10  for 
the  mildly  deformed  shapes  of  Figures  7 and  11  to  around  30  for  the  deeper  and  re-entrant 
grooves  in  Figure  13.  The  orthogonality  weight  for  the  smoothing  functional,  W2j  was  kept 
at  10  for  all  grids.  The  Jacobian  weight,  ^^;l,  ranged  between  .01  and  .2  for  the  grids  in 
Figure  7,  but  was  set  to  .5  for  the  deeper  deformations  shown  in  Figures  11b,  12,  and  13. 
The  grids  were  able  to  maintain  a significant  amount  of  orthogonality  and  smoothness  both 
within  the  interior  and  along  the  boundary  as  the  grid  points  redistributed  themselves  to 
follow  the  interface. 


5.  Conclusions 

The  development  of  an  algebraic  grid  generation  system  to  track  a solid-liquid  interface 
has  been  discussed.  The  proposed  grid  generation  mapping,  composed  of  tensor  product 
B-splines,  effectively  approximates  interface  shapes  of  varying  degrees  of  complexity.  For 
each  shape  the  initial  coefficients  are  chosen  to  approximate  a transfinite  blending  function 
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Figure  11:  Grids  show  the  mildly  deformed  sinusoidal  shapes  commonly  observed  during  the 
initial  stages  of  the  deformation  of  a solid-liquid  interface. 


Figure  12:  Grids  show  deeper  grooves,  but  orthogonality  is  mcuntained  at  the  interface. 
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(a) 


(b) 


Figure  13:  Grid  on  left  shows  deep  and  narrow  groove.  Grid  on  right  shows  deep  re-entrant 
bulb-like  cell. 

mapping  that  produces  nearly  orthogonal  grid  lines  near  the  interface.  Adjusting  the  coef- 
ficients to  minimize  a discrete  smoothing  functional  smoothes  the  spacing  of  the  grid  lines, 
untangles  overlapping  lines,  and  enhances  the  orthogonality  of  the  initial  grid  even  as  the 
interface  becomes  more  distorted.  Smoothness  and  orthogonality  is  improved  along  the  outer 
boundary  by  adjusting  those  boundary  coefficients  whose  movement  will  not  alter  the  shape 
of  the  boundary. 

Currently,  the  grid  generation  system  allows  some  interaction  as  the  interface  changes  so 
that  orthogonality  and  smoothness  parameters  can  be  adjusted  as  the  interface  cell  deepens. 
Ideally,  this  should  all  be  done  automatically  by  the  code.  Additional  study  will  be  continued 
in  this  area.  The  next  phase,  which  has  already  begun,  is  to  couple  the  grid  generation 
algorithm  with  equations  that  deterrmne  the  interface  shape  so  that  the  system  can  be  used 
in  the  numerical  analysis  of  microstructures  that  develop  during  directional  solidification. 
Although  the  system  will  be  designed  to  track  the  deformation  of  a planar  interface  into  a 
deep  ceU,  it  may  also  be  used  to  improve  calculations  in  phase  field  models  [20,  21,  22]  where 
the  interface  is  not  viewed  as  a curve  or  surface  with  zero  thickness.  In  such  models  accuracy 
can  be  improved  by  concentrating  the  grid  points  in  the  area  neaj  the  interface  even  if  the 
interface  is  not  tracked  exactly. 

Disclaimer 

Identification  of  commercial  products  in  this  paper  does  not  imply  recommendation  or  en- 
dorsement by  NIST. 
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